Simulation of particle mixing in turbulent channel flow due to intrinsic fluid velocity 

fluctuation 
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We combine a DEM simulation with a stochastic process to model the movement of spherical 
particles in a turbulent channel flow. With this model we investigate the mixing properties of two 
species of particles flowing through the channel. We find a linear increase of the mixing zone with 
the length of the pipe. Flows at different Reynolds number are studied. Below a critical Reynolds 
number at the Taylor microscale of around R c ~ 300 the mixing rate is strongly dependent on the 
Reynolds number. Above R c the mixing rate stays nearly constant. 

PACS numbers: 47.27.T-, 45.70.Mg, 83.50.Xa, 47.27.E- 



I. INTRODUCTION 

Transport of granular material by a fluid plays an im- 
portant role in many environmental and technological 
processes. In nature examples are sand storms in the 
desert, the formation of dunes [l], and river deltas 
3 or the dispersion of particles like pollen or carbon 
black in the atmosphere or microorganisms or pollution 
in the ocean 0] or also particle transport and separa- 
tion in lung- like structures [H and rough channels Q- In- 
dustrial and technical applications range from transport 
of contaminant and pollution, to filtration Q, combus- 
tion, deposition of sprays, jet mills, fluidized beds, pneu- 
matic transport of powders, etc. Together with trans- 
port of particles naturally mixing and segregation oc- 
curs. In many industrial applications a fine control over 
these two effects is essential @. To achieve this, differ- 
ent procedures have been developed. Mechanical mixers 
like tumbling mills, drum mixers or impeller type mixer 
are widely used. Another successfully applied technique 
consists of using jet nozzles that cause bubble formation 
in fluidized beds. In this situation the mixing of particles 
is driven by inhomogeneitics in the particle density. 

In many industrial applications pneumatic conveyor 
systems are used to transport particles from one pro- 
cess section to the next. Such conveyor systems usu- 
ally consist of a channel wherein a turbulent air flow is 
used to fluidize and transport particles. In these sys- 
tems the particle densities can be considered to be ho- 
mogeneous enough to avoid mixing driven by density in- 
homogencities as mentioned above. Still mixing can be 
observed and this effect could be used in industrial appli- 
cations in the future. In this situation the dispersion of 
particles is determined entirely by the fluid flow. In flows 
with a high Reynolds number the statistics of velocities 
and accelerations present inside the turbulent flow are 
well known from experiments 0] . We study the influence 
of these intrinsic fluctuations on the mixing of systems 



of constant particle density. For this we couple a parti- 
cle simulation to a stochastic model for the fluid velocity 
fluctuations [l(| that reproduces well the observed statis- 
tics of the fluid velocity and acceleration. We investigate 
the influence of the strength of the turbulence on the 
mixing. 

Different approaches for numerical simulations of dif- 
ferent kind of mixers have been developed (see e.g. 
Bcrtrand et al. 
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for a review) . A popular type of mod- 
els simulates the particles individually and treats occur- 
ring collisions explicitly (r34T^ | . For simulating fluidized 
beds one usually couples a particle based method like the 
Discrete Element Method (DEM) Q with a fluid model. 
A good review of current techniques has been published 
by Deen et al. [l5| . The basic ideas for simulating bub- 
ble formation were first introduced by Tsuji et al. fl(| 
and Hoomans et al. (l7j or was more recently also used 
in [H]. 

A large variety of models for fluid-particle coupling 
have been proposed and a good overview of their cur- 
rent state can be found in the review of Zhu et al. fl9j . 
One important class of techni ques consists in using direct 
numerical simulation (DNS) [20( or Lattice-Boltzmann 
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(LB) [21| to solve the Navier-Stokes equations and cou- 
ple the resulting fluid velocity field to the particle motion 
by an empirical drag law j22[. Unfortunately these tech- 
niques are limited to only moderate values of Reynolds 
numbers. For higher Reynolds numbers one has to re- 
lay on turbulence models like Large-Eddy simulations, 
mixing length models or other turbulence models. These 
models are usually not able to resolve the smallest fluctu- 
ations of the fluid field but rather calculate a (temporal 
or spatial) average of the fluid velocity field [23[ . 

In this paper we investigate the influence of the intrin- 
sic fluctuations of the turbulent velocity field on the mix- 
ing. We study two species of spherical particles flowing 
through a rectangular pipe as sketched in Fig. [1] We ne- 
glect all influences from the wall on the fluid velocity field 
but instead use a stochastic model [l(| that reproduces 
well the measured acceleration and the velocity statistics 
of tracer particles in fully developed turbulence by La 
Porta et al. 0. We start first with a two dimensional 
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FIG. 1. (color online) Schematic side view of a rectangular 
channel of width W and height H. Two species of particles 
are transported through the channel by a turbulent flow of 
mean velocity (u s ). The simulation concentrates on a small 
slice of size L x x L y located at the center of the channel. The 
two species enter the simulation on the left side and leave it 
again on the right side. During this displacement the particles 
are mixed due to the turbulent flow. The upper part of the 
figure shows a snap-shot of the system seen from above for 
i?A = 1000 and p v = 0.5. 
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acts on particle i. Here A; is a spring stiffness, 70 a viscous 
damping coefficient and the relative velocity of the 
particles in the normal direction 



(6) 



B. Calculation of the drag forces 



Next we have to specify how the drag force is calcu- 
lated. Given a (local) fluid velocity u s we follow Bini and 
Jones [25I ] by assuming the simplified expression for the 
acceleration of particle i due to the flow 



U s - Vi 



a, 



(7) 



system and extend it to three dimensions. The influence 
of different Reynolds numbers as well as different densi- 
ties is analyzed. 



II. THE MODEL 

Our model is based on a two-dimensional DEM sim- 
ulation (2~H in a rectangular domain. The particles are 
non-dcformable disks with radius and density pi. The 
equations of motion for a particle i are given by 



(1) 



where rrii and Xj are the particle mass and position and 
Fj is the sum of all forces acting on the particle. In our 

(c) 

model Fj is the sum of a collision F,- and a drag force 



r(d) 



due to the fluid in the channel, 



,(c) 



r (d) 
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Given Fj, the equations of motion are integrated using 
the standard Verlet algorithm. 



where Vj is the particle velocity and t is a relaxation 
time. In a laminar flow one can choose r to be constant 



resulting in a drag force 



.(d) 



acting on particle i pro- 



portional to the velocity difference 
r (d) 



cx u s 



(8) 



In the turbulent case one generally assumes a quadratic 
dependence of the drag force on the velocity difference 



F„ (d) oclus-v, 



(9) 



To achieve this an appropriate expression [251 ] for r is 
_ 1 _ 3 p/C D . 



8 pm 



(10) 



where pf, pi, ri and Cd are the fluid density, the particle 
density, the particle radius and the drag coefficient, re- 
spectively. For spherical particles at high Reynolds num- 
bers C D « 0.45. 

At this point it is important to note that the fluid is not 
influenced by the particles. In a more detailed simulation 
one would need to include a feedback on the fluid due to 
the energy transfer from the fluid to the particles. 



A. Calculation of the collision force 

To model the collision force between two particles with 
positions Xj and x 3 , we first calculate their overlap 



Sij = (n+rj) - |xj 



(3) 



C. Determination of the turbulent flow 

To complete our model we need to define the turbulent 
flow velocity u s . We first split u s into the mean stream 
velocity (u s ) and a stochastic part u t : 



If Sij is positive, a repulsive force 



u t . 



(11) 
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Calculating the mean flow velocity (u s ) is rather easy 
and would also account for the wall effects. For example 
one could apply a turbulent viscosity model [23| , like the 
mixing length model or the widely used k — e model, to 
determine (u s ). 

Determining the stochastic part u t of the flow is more 
difficult. Several measurements [t| H|| [27j have shown a 
highly non-Gaussian behavior of the Lagrangian acceler- 
ation distribution of fluid particles. Sawford [28[ intro- 
duced a system of stochastic differential equations (SDE) 
for the time evolution of the velocity and acceleration 
of a fluid particle. This model correctly reproduces the 
observed Gaussian distribution of the velocities. Unfor- 
tunately it fails to show the highly non-Gaussian distri- 
butions for the acceleration. To overcome this problem 
Pope and Chen [2!| introduced an Uhlenbeck-Ornstein 
process for the evolution of the (now fluctuating) energy 
dissipation rate e. Reynolds [10| combined these ideas 
and formulated a system of SDEs for the logarithm of the 
dissipation rate x = l°g( £ / ( £ ))j an d the components of 
the accelerations at and velocities u t of tracer particles in 
fully developed turbulence. This model reproduces well 
[30l ] the observed probability density functions as well as 
the autocorrelation functions of the velocity and acceler- 
ation: 



dx 
dat 



(x - (x^T-'dt + y / 2^p7 I ^ 1 (12a) 



= - TV 1 + t,- 1 - g 



_i d(T a t \e 
a *l e dt 



atdt 



-T L t„ 1 utdt 



2<jl(T- 1 +t^)T-% 1 d^ 



du t = a t dt. 
In Eq. (|12a|) the variance of 



(12b) 
(12c) 



a l x = -0.354 + 0.289 log R x 
given by (x) = 

2 



is approximated by 
its mean value is 
0.5c 2 , the relaxation time scale 



T x = 2ct„/(Co (e)) and (e) is the mean energy dissi- 
pation rate. The Reynolds number at the Taylor mi- 
cr oscalc R\ is related to the Reynolds number by R\ — 
a/15Rc. In Eq. (|12b[) we have the energy-containing 
scale Tl = 2<7 2 /(Coe), the energy-dissipation scale t v = 
Cqv 1 ! 2 /(2ao£ 1 ^ 2 ), the conditional acceleration variance 
a a t \s = ao^ 3 ^ 2 ^ 1 ^ 2 1 two universal Lagrangian velocity 
structure constants clq = 3.3 and Co = 7.0, the kinematic 
viscosity v and the velocity variance <r 2 . The values of ao 
and Co are determined by demanding consistency with 
Kolmogorov's (1941) hypothesis (32| and fitting them to 
experimental data |10(. Finally there are two indepen- 
dent Wiener processes, i.e. Gaussian distributed random 
numbers, d£i and d^2 with zero mean and variance dt. 

With this model we are able to calculate each com- 
ponent of the turbulent velocity u t . We combine one 
equation for x and two equations for at and u t to get a 
two-dimensional vector u t for the velocity of a tracer par- 
ticle. We then generate a set of these tracer particles and 
evolve them in time according to Eq. 1121 When a new 



particle enters the system, we randomly chose a tracer 
particle out of this set and use this (tracer) particle to 
calculate the stochastic forces acting on the (real) parti- 
cle until the particle leaves the system. This means that 
we add the mean velocity (u s ) and the stochastic part u t 
from the assigned tracer particle to get a "stream line" 

which we then use to calculate the drag force F 



(d) 



III. SIMULATION SET-UP 

We applied the model described in the last section to 
a system with two types of particles flowing through a 
channel. A sketch of the system is shown in Fig. [TJ Both 
types of particles have the same radius r = 10 -5 m and 
density p = 3.0 • 10 3 kg m -3 (but a different "color" 
to distinguish them). The spring stiffness and viscous 
damping coefficient are set to k = 1.0 kgs -2 and 70 = 
1.0 • 10~ 6 kgs -1 and kept constant for all simulations. 
The fluid inside the channel has a kinematic viscosity 
v = 1.64 • 10~ 5 m 2 s _1 . The channel has a width W 
and height H. In our case we have W = 0.2 m and 
H = 0.358 m. To get a two dimensional system we first 
cut out a horizontal slice from the channel. We further 
want to simplify the setup and study the basic effects 
of the stochastic part in Eq. (fTTj) . Therefore we con- 
centrate on a small region in the center of the channel 
with dimensions L x = 0.025 m and L y = 0.002 m. Next 
we set the mean velocity to be constant throughout the 
whole system and pointing into the positive x-direction, 
(u s ) = ((u s ) ,0). This of course ignores wall effects, be- 
cause the mean flow profile in a channel is not constant 
and drops to zero at the boundaries. But as mentioned 
before our simulation focuses only on the region in the 
center of a large channel where the mean velocity is ap- 
proximately constant. These assumptions influence the 
choice of our boundary conditions which will be discussed 
later. 

The magnitude of (u s ), as well as the mean energy 
dissipation rate (e) and the velocity fluctuations g\ de- 
pend on the Reynolds number. The values used in our 
calculations are listed in Tab. U These parameters were 
determined in the following way: 

1. For pipe flow the Reynolds number at the Taylor 
microscale [23| can be calculated by 



Rx = \/15^£> H . 



(13) 



Here Du is the hydraulic diameter. For a rectan- 
gular duct of width W and height H this is given 

by 



D H = 



2WH 
W + H' 



(14) 



With this expression one can calculate (u s ) for a 
given Reynolds number R\. 
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TABLE I. List of Reynolds number dependent parameters in 
our simulations. 
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2. In fully developed turbulence there is an alternative 
expression for the Reynolds number [23j : 



R\ = 




(15) 



where A = ^JVovoWj (e) is the Taylor microscale. 
For a fixed Reynolds number R\ this gives us a re- 
lation between the mean energy dissipation rate (e) 
and the velocity fluctuations a u . Another relation 
between these two quantities is given by 



L ' 



(16) 



where L is a "typical" length scale. We assume that 
L is of the same order as £>h, i-e. L ~ Z?h- This 
allows us to calculate values for (e) and a u . 

3. The Stokes number, quantifying the effect of inertia 
on the suspended particles, can be calculated [H| 
as 



St = 



n 



where r, is the particle relaxation time 



18vpf 

and is the Kolmogorov time scale 



(17) 



(18) 



(19) 



Next we have to specify the boundary conditions of 
our system. New particles enter the system at the left 
side and leave it again at the right side (See Fig. [1]). 
The boundary conditions in the transverse direction are 
more complex: Particles are elastically bounced back by 



reflecting their velocity vector such that the vertical com- 
ponent points back into the system. This leaves the total 
kinetic and potential energy of the particles constant. 
Additionally, we exchange the to the particle assigned 
tracer particle that is used to calculate the stochastic 
part of the fluid velocity. This is done by drawing uni- 
formly distributed a tracer particle out of our set of all 
tracer particles. Such a boundary condition is reason- 
able because we concentrate only on a small region in 
the middle of the channel and our boundaries are not the 
walls of this channel. One can imagine that whenever 
a particle hits our boundaries, in reality it would leave 
our domain of simulation. So, to keep the (local) particle 
number constant, we think of the reflected particle as a 
new particle that would enter our simulation domain at 
the same place and time where the old particle left. 

Finally we have to specify the density of our system by 
fixing the volume fraction p v that is occupied by particles 



Pv 



(ni ■ V\ + ri2 ■ V 2 ) 
L x ■ L,. 



(20) 



Here rik and Vk denote the number of particles of type k 
and their respective volume. This determines how many 
particles are on average in the system. Using the mean 
stream velocity we can estimate how many particles leave 
the system and therefore how many new particles we have 
to insert at the opposite side into the simulation to keep 
p v constant. We divide the left boundary of the system in 
two separate equally sized intervals, insert the one species 
into the right and the other species into the left interval. 
New particles are placed in such a way that they do not 
overlap with existing ones. A section of the system at the 
end of the channel is shown in the upper part of Fig. [T] 



IV. RESULTS 

To quantify the quality of mixing, we cut out vertical 
slices of a certain width from the system. These slices we 
divide into horizontal intervals and count how many par- 
ticles of which type are in each interval. Averaging over 
many time-steps we get a time-averaged relative density 
for each interval. To be more specific, let's look at a slice 
located at x s . We divide this slice into n s horizontal 
intervals. This gives us n s small boxes Bi, centered at 

(x s , ^ ■ (i + 0.5)^. During the simulation we count the 
number of particles of each type, namely TV^.i particles of 
type 1 and particles of type 2 in the box From 
this we can calculate the relative densities 



N, 



Pi,k 



k = 1.2. 



(21) 



The behavior of these relative densities for the two 
species of particles in the center of the channel along one 
slice is shown in the inset of Fig. [21 
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A. Growth of the mixing zone 

Now we will quantify the size of the mixing zone. 
Therefore we calculate 

p% = i - \Pi,i - Pi,2\ ( 22 ) 

to get a quantity which is close to zero in regions where 
the mixing is bad, increases the more the particles are 
mixed and approaches unity when the two species are 
well mixed. Examples for the behavior of //, along two 
slices arc shown in Fig. [2] 

We now calculate the area A below the measured points 
fii to get a quantity that specifies the size of the mixing 
zone 

n s 

A = X>- (23) 
i=i 

We can calculate A for every slice and find a lin- 
ear increase with increasing length of the channel (see 
Fig. [3] (a)). By altering the Reynolds number we find an 
increase of the slope of A(x s ) with increasing Reynolds 
number. Above a certain Reynolds number R c the qual- 
ity of mixing does not increase any more, i.e. for R\ — > oo 
the slope dA/dx approaches a constant value. This re- 
sult is shown in Fig. 01 In our simulation R c sa 300. This 
result has implications on industrial applications where 
mixing in turbulent flows is of key importance. 

If we fix the Reynolds number and alter the density, we 
find a faster mixing for lower densities (see Fig. [3] (b)). 
The dependence of the growth rate of the mixing zone on 
the density is shown in the inset of Fig. [3] (b). With our 
current implementation it is possible to reach densities 
of about 0.6, because particles are not allowed to overlap 
when they enter the system. By using another way of 




FIG. 2. (color online) Example of the behavior of fii along 
two slices with R\ = 1000 and p v = 0.2. The inset shows 
the behavior of the relative densities of the two species in the 
center of the channel at x„ = 0.007625 with Rx = 1000 and 
Pv = 0.2. 




position inside the channel x 



FIG. 3. (color online) Behavior of the size A of the mixing 
zone in two dimensions for different Reynolds numbers and 
densities: (a) Behavior of A(x) for constant density p v = 0.2 
and different Reynolds numbers R\. One observes a linear 
increase of the mixing zone along the channel. The mixing 
gets faster with increasing Reynolds numbers until it reaches 
a limit, i.e. the slopes dA/dx of the linear fits stay constant. 
Curves are shifted vertically to better distinguish the indi- 
vidual curves, (b) Behavior of A(x) for constant Reynolds 
number R\ — 800 and different densities. Systems with lower 
densities exhibit faster mixing. The inset shows the growth 
rate dA/dp of the mixing zone. 



inserting particles into the system it may be possible to 
reach even higher densities close to the value for random 
close packing p vc w 0.84. Because particles can be elas- 
tically deformed, we expect that even at densities close 
to random close packing we can still observe mixing and 
the growth rate of the mixing zone goes to a finite value. 
This tendency can be seen in the inset of Fig. [3] (b). 



B. Lateral dispersion 

As an alternative to measuring the size of the mixing 
zone we use a dispersion relation to quantify the rate of 
mixing. The motion of tracer particles can be character- 
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FIG. 4. (color online) Dependence of the slope dA/dx on 
the Reynolds number R\. Above a certain Reynolds number 
R c ^ 300 the slope stays almost constant. Empty symbols 
represent the results in 2D and filled symbols are the ones for 
3D, discussed in section IVl 



ized [34| by a convection-diffusion equation 

dpk{ * ,£) +u s (x)-Vp fc (x,t) = V- (D(p fc ,x)Vp fc (x,t)), 

(24) 

where pk (x, t) is the relative density of particles of 
species k at position x and time t, u s (x) the fluid ve- 
locity and Z?(pfe,x) the coefficient of dispersion. Eg. |2"41 
can be simplified as follows: The first term vanishes by 
assuming a steady state 



at 



= o. 



(25) 



The second term simplifies by replacing the fluid velocity 
by the mean fluid velocity 



u s (x) = (u s ) =((«b>,0) 
The right hand side simplifies by assuming 

D(pk,x) = D = const, 
and by using according to Ref. [H[ 



d 2 p k ( x ,t) <<; d 2 Pk {^t) 



Therefore Eq. 



9a; 2 
| simplifies to 



dy 2 



(u s ) 



dp k „d 2 p k 



dx 



D- 



dy 2 



(26) 
(27) 

(28) 
(29) 



If one assumes Heaviside functions as the initial density 
profiles at x = 0, the solution of Eq. [29] is given [3f| by 



ph = \ l±erf _JL 




(30) 



where erf (x) is the Error function. Rewriting this equa- 
tion leads to 



erf" 1 (1 - 2 • p fc ) = 



y 



2VD yjx/ (li s ) 



(31) 



For every position of a slice x = x s we can now calculate 
the left and right hand side of Eq.[2Ualong the y-direction 
and fit for the dispersion coefficient D. In theory this 
should lead to the same dispersion coefficient for all slices. 
As an example we show in Fig. [5] the results for a system 
with p v = 0.2. As can be seen D is not constant along 
the channel and depends on the position x inside the 
channel and the Reynolds number R\. Rescaling the 
values of D(x, R\) by 1/D(x — 0.02, R\) leads to a data 
collapse of the curves. The inset of Fig. [5] shows that 
the values of D(0.02,i?A) seem to depend quadratically 
on the Reynolds number R\. The same behavior can 
be observed for other densities as well. In Fig. [5] one 
may also suspect a crossover between two linear regimes. 
However, this does not seem to be observable for other 
densities. 

For small values of x the variation of D may be due 
to the initial conditions of particles entering the system, 
where all particles have the same velocity (u s ) pointing in 
the positive a;-direction. But because the values of D do 
not converge towards a constant value it seems likely that 
the lateral dispersion of particles in our system cannot be 
described by normal dispersion. To conclusively judge 
whether our system shows anomalous dispersion or not 
further investigations e.g. of the behavior of the mean 
squared displacement of one particle over time would be 
needed. 

Even the dispersion coefficients being not constant 
throughout the system, we tried to compare them with 



Q 




0.02 



FIG. 5. (color online) Data collapse of the dispersion coeffi- 
cients D(x, Rx) for p v = 0.2. For every Reynolds number R\ 
the corresponding curve along the channel was rescaled with 
the dispersion coefficient at x — 0.02. The inset shows the 
quadratic dependence of the values D(x = 0.02, R\) on the 
Reynolds number R\. 
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the turbulent diffusivity of Taylor [34j, |36[ defined by 

D E = alT L . (32) 

Inserting the formula for the Lagrangian integral time 
scale Tl given in Sec. Ill Cl leads to 



D E = 



2u 
15Co" 



Rl 



(33) 



This explains the observed quadratic dependence of D 
on R\ . Inserting the values for the viscosity v and model 
parameter Co results in values for De of the order 1CP 4 
for small Reynolds numbers to 10 -1 for large ones. These 
values are clearly orders of magnitude above the values 
observed in our simulations. 



C. Particle-particle and fluid-particle interactions 

The motion of particles in our system is determined 
by two types of interactions: The particle-particle colli- 
sions and the drag force from the fluid. The initial con- 
ditions of particles entering our system are chosen such 
that these particles do not overlap and their velocities 
are equal to the mean flow velocity (u s ). If we remove 
the stochastic part of u s particles would just fly along 
straight lines and one would not observe any mixing at 
all. By adding random components to the initial velocity 
of the particles the mixing would be similar to that in a 
granular gas. In our simulation the particles all have the 
same velocity when entering the system and variations of 
the initial velocity are introduced by the stochastic forces 
described in Sec. Ill B land III Cl To estimate the influence 
of the particle-particle and fluid-particle interactions we 
can compare the mean free time between collisions in a 
granular gas [37| 



2y/2 Pv 



(34) 



with the Kolmogorov time of the turbulent part u t of 
the fluid velocity. The mean free time r g is related to 
the mean free path A of a granular gas by dividing A by 
the mean particle velocity. Inserting the parameters for 
our simulations shows that the Kolmogorov time is about 
10 — 1000 times larger than r g . This means that many 
particle collisions occur between changes of the fluid ve- 
locity. In this case the Stokes number plays an important 
role. Tab. U shows that for R x < 400 the Stokes num- 
ber is smaller than one which means that particles fol- 
low the fluctuations in the fluid velocity u s rather closely 
and the particle collisions can be seen as small pertur- 
bations only. Therefore the actual values of the spring 
stiffness k and viscous damping coefficient 70 should not 
have a significant influence on our results in this range 
of Reynolds numbers, even for the limiting case of hard 
spheres (k — > 00). However, St > 1 for R\ > 400 which 
means that particles will be increasingly less influenced 



by changes in the fluid velocity and particle-particle col- 
lisions become more important. This may be a reason 
for the saturation of the mixing rate for large Reynolds 
numbers observed in Sec. IIV Al Considering variations 
of k and 70, in particular the case of hard spheres, may 
have an influence on this asymptotic value. 

It is also worth noting that numerous publications |3S 
|4~0 | report that particles in turbulent flows tend to con- 
centrate in regions of low vorticity and high strain rate. 
Bee et al. showed that this effect of preferential con- 
centration can be observed for heavy particles and Stokes 
numbers St < 1. From Tab. U one can see that this crite- 
rion is fulfilled for R\ < 400 and therefore clustering of 
particles should be observable in our simulations. Indeed 
when looking at a system with zero mean flow velocity 
(u s ) = and periodic boundary conditions clustering of 
particles was observed. In the systems studied in this 
paper one may recognize such clustering at the end of 
the channel. The reason for this less obvious occurrence 
of preferential concentration is probably that the initially 
homogeneous distributed particles leave the system again 
before a clear clustering can be observed. 



V. EXTENSION TO THREE DIMENSIONS 

We also extended our model to three dimensions. Most 
parts of the two-dimensional model carry over to the 
three-dimensional one in a straight forward way. Now, 
particles become spheres instead of disks. The collision 
and the drag force are calculated exactly the same way 
by just replacing two-dimensional vectors by their cor- 
responding three-dimensional ones. The stochastic part 
u t of the flow velocity u s is calculated by adding a third 
component of a t and Ut to the computation. Finally we 
extend our formerly two-dimensional domain of simula- 
tion in the z-direction by a length L z to get a three- 
dimensional rectangular channel located in the center of 
the big channel. Particles still enter the system at the left 
side and leave it at the right side. The boundary condi- 
tions on the other four sides of the channel are the same 
as in the two-dimensional model: Particles are bounced 
back from the wall and their attached tracer particles 
are exchanged. We studied a system of size L x = 0.01 m, 
L y = 0.001 m, L z = 0.0003 m and two different densi- 
ties p v = 0.2 and p v = 0.5. This system is smaller than 
the two-dimensional one because the computational cost 
for a three dimensional computation is much higher. All 
other parameters stay the same as in two dimensions. 

To quantify the quality of mixing we used a simi- 
lar approach as in the two-dimensional case. Again 
we cut out slices perpendicular to the x-axis at differ- 



ent positions x s . Every slice we divided into n 
tervals along the y-axis and n s . z intervals along the z 
axis. This gives us n s , y ■ n s<z boxes Bij centered at 

(x s , + 0.5), ^7 • (j + 0.5)) . During the simula- 

tion in every box Bij we again counted the number of 
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FIG. 6. (color online) Behavior of the size A of the mix- 
ing zone in three dimensions for different Reynolds numbers 
at p v = 0.2. One again observes a linear increase of the 
mixing zone along the channel. The mixing gets stronger 
with increasing Reynolds numbers until it reaches a limit, i.e. 
the slopes dA/dx of the linear fits stay constant. Curves are 
shifted vertically to better distinguish the individual curves. 



particles of every species, i.e. Nij i and Nij %. Because 
the system has a translational symmetry along the z-axis 
we calculated 



(Ni, 



i,j,ki 



1,2. 



(35) 



This gives us the mean particle number of every species 
along the y-axis. From this we again can calculate the 
relative densities 



PiM — 



i,k 



(N iA ) + (N ii2 ) 



k = 1,2. 



(36) 



The qualitative behavior of these densities is the same 
as in the two-dimensional model shown in the inset of 
Fig. [2] From here on the calculation is exactly the same 
as before. We first calculate fa, then compute the area A 
below these curves and plot A versus the position x in- 
side the channel. This plot is shown in Fig. [5] for p = 0.2 
at different Reynolds numbers. Again we find a linear 
increase with respect to the length of the channel. By 
changing the Reynolds number we see a similar increase 
of the slope of A(:r) as in two dimensions. To compare 
the two and three dimensional models we plot in Fig. [4] 
the slopes of A(x) for different Reynolds numbers R\ and 
different densities p v for the two- and three-dimensional 
cases. As expected the rate at which the mixing zone 
grows increases with higher Reynolds numbers and goes 



to a constant value when R\ — > oo. The mixing is also 
slower for higher densities. More interesting is the com- 
parison between the mixing rate in two and three dimen- 
sions for a fixed density: One would expect, that the mix- 
ing is faster in three dimensions because particles have 
an additional degree of freedom in three dimensions and 
can move over each other more easily. This behavior we 
can observe for high Reynolds numbers. But surprisingly 
for small Reynolds numbers the mixing in three dimen- 
sions is actually slower than in two dimensions. Another 
interesting result is that the nearly constant value of the 
mixing rate is approached faster in two dimensions than 
in three. 



VI. CONCLUSION 

By combining a molecular dynamics simulation with a 
PDF method for modeling turbulent flows, we were able 
to develop a simple model to investigate the mixing of 
two species of particles inside a channel. The mixing is 
caused solely by the stochastic force measured in fully 
developed turbulence. We found that the mixing zone 
increases linearly with the length of the channel. The 
mixing is stronger at higher Reynolds numbers and at 
lower densities. For small Reynolds number the mixing 
is stronger in two dimensions than in three. We further 
investigated the lateral dispersion of particles and found 
hints toward anomalous dispersion. 

Further improvements of our model would include a 
more detailed modeling of the mean fluid flow in the 
channel to include the wall effects. At the same time the 
effect of gravity should be included in the simulations. 
One also should introduce a feedback from the particles 
back to the fluid and also different sizes and shapes of 
particles could be studied. Here we investigated non- 
dcformablc particles. In future work it would be inter- 
esting to study the influence of particle deformations on 
the mixing. Whether or not anomalous lateral dispersion 
can be observed remains an open question. Finally in- 
vestigating situations where the mean free time between 
particle collisions becomes of the same size as the turbu- 
lent velocity autocorrelation time could also significantly 
change the mixing behavior. 
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